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The full kinetic dynamics of a perpendicular collisionless shock is studied by means of a one-dimensional 
electromagnetic full particle simulation. The present simulation domain is taken in the shock rest frame in contrast 
to the previous full particle simulations of shocks. Preliminary results show that the downstream state falls into a 
unique cyclic reformation state for a given set of upstream parameters through the self-consistent kinetic processes. 
Key words: Collisionless shock; particle-in-cell simulation; shock rest frame 



1. Introduction 

Collisionless shocks are universal processes in space and 
are observed in laboratory, astrophysical, and space plasmas, 
including astrophysical jets, an interstellar medium, the he- 
liosphere, and the planetary magnetosphere. The classic pic- 
ture of collisionless shocks is that they arise in a collision- 
less plasma streaming with a supersonic velocity against an- 
other collisionless plasma. The generation of plasma waves 
and turbulences, plasma heating and acceleration, and elec- 
tromagnetic radiation processes at collisionless shocks are 
of great interest in astrophysics, solar physics, planetary 
physics, geophysics, and plasma physics. Detailed scien- 
tific issues in collisionless shock physics are described in 
Lembege et al. (2004). Kinetic simulations of collisionless 
shocks are essential approaches to solve these scientific is- 
sues. 

There are several different methods for exciting collision- 
less shocks in numerical simulations of plasmas. The mag- 
netic piston method (Lembege and Dawson, 1987a; Lem- 
bege and Savoini, 1992) is widely used in full particle simu- 
lations, in which a plasma is accelerated by an external cur- 
rent pulse applied at one side of the simulation domain. The 
simulation domain is taken in the upstream rest frame. The 
plasma is pushed by the "magnetic piston" into the back- 
ground plasma, and the external pulse develops into a shock 
wave. Another method widely used is the injection method 
(Quest, 1985; Shimada and Hoshino, 2000; Schmitz et al, 
2002a, b; Nishimura et al, 2003; Scholer et al, 2003; Lee 
et al, 2004), in which a plasma is injected from one side of 
the simulation domain and is reflected back when it reaches 
the other side. (Therefore this method is also called the re- 
flection method or wall method.) The simulation domain is 
taken in the downstream rest frame, and a shock wave is ex- 
cited due to the interaction between the reflected and the in- 
jected plasma. Additional methods, such as the flow-flow 
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method (Omidi and Winske, 1992) and the plasma release 
method (Ohsawa, 1985), are used in hybrid and full particle 
simulations. A detailed description of each method is given 
in Lembege (2003). 

In these methods, an excited shock wave propagates up- 
stream. Therefore, it is necessary to take a very long simula- 
tion domain in the propagation direction of the shock wave 
in order to study a long-time evolution of the excited shock 
wave. This makes it difficult to perform multidimensional 
simulations even with a present-day supercomputers. In this 
letter, we first attempt to perform a full particle simulation of 
a collisionless shock in the shock rest frame: the collision- 
less shock is excited by using the "relaxation method" which 
was used in hybrid code simulations in 1980s (Leroy et al, 
1981, 1982; Kan and Swift, 1983). This method has not been 
used in full particle simulations due to several difficulties in 
numerical techniques. 

2. Simulation Model 

Our simulation code (Umeda, 2004) is an improved 
version of one-dimensional Kyoto university electromag- 
netic particle code (Omura and Matsumoto, 1993), where 
Maxwell's equations and the relativistic equations of mo- 
tion for individual electrons and ions are solved in a self- 
consistent manner. The continuity equation for charge is also 
solved to compute the exact current density given by the mo- 
tion of the charged particles (Umeda et al., 2003). 

The simulation domain is taken in a one-dimensional sys- 
tem along the s-axis. The initial state consists of two uni- 
form regions separated by a discontinuity. In the upstream 
region that is taken in the left hand side of the simulation 
domain, electrons and ions are distributed uniformly in x 
and are given random velocities (v x , v y ,v z ) to approximate 
shifted Maxwellian momentum distributions with the drift 
velocity u x \, number density n\ = eom e ujp el /e 2 , isotropic 
temperatures T e \ = m e v^ el and Tn = rriiV^, where 
m, e, ojpi, and vti are the mass, charge, upstream plasma 
frequency, and upstream thermal velocity, respectively. In 
this letter, subscripts "1" and "2" denote "upstream" and 
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"downstream", respectively. The upstream magnetic field 
Boyi = —m e ui ce i/e is also assumed to be uniform, where 
u> c i is the upstream cyclotron frequency (with sign included). 
The downstream region taken in the right-hand side of the 
simulation domain is prepared similarly with the drift veloc- 
ity u x2 , density n 2 , isotropic temperatures T e2 and T i2 , and 
magnetic field Bq V 2- In this letter we assume a perpendicular 
shock (i.e., Bq x = 0). As a motional electric field, an exter- 
nal electric field Eq z = —u x iBo v i — —u x iB§ y i is also 
applied in both the upstream and downstream regions. At 
the left boundary of the simulation domain, we inject plas- 
mas with the same quantities as those in the upstream region, 
while plasmas with the same quantities as those in the down- 
stream region are also injected from the right boundary. We 
adopted absorbing boundaries to suppress the non-physical 
reflection of electromagnetic waves at both ends of the sim- 
ulation domain (Umeda et al, 2001). 

In the present simulation, the time, velocity, and posi- 
tion are normalized by the initial upstream electron plasma 
frequency Lo pe \ = 1.0, upstream electron thermal velocity 
vtei = 1-0, and upstream electron Debye length A e i = 
vtei/^pei = 1.0, respectively. The initial temperatures in 
both upstream and downstream regions are assumed to be 
isotropic. In the upstream region, we assume a low beta and 
weakly magnetized plasma, such that j3 e x = fin = 0.125, 
and u! ce i/uj pe i ~ —0.05, which are similar to the recent full 
particle simulations (Shimada and Hoshino, 2000; Schmitz 
et al, 2002a, b; Lee et al, 2004). The light speed in the 
present simulation is given as c = 80.0. The bulk flow ve- 
locity of the upstream plasma is assumed to be u x \ — 4.0, 
which corresponds to the Alfven Mach number Ma = 10.0. 
The ion-to-electron temperature ratio, rx = Ti/T e is as- 
sumed to be unity in the upstream region (r^i = 1). In 
addition to the upstream quantities u x i, ^ P ei, ^cei, Vtei, 
and rxi, we need the downstream ion-to-electron tempera- 
ture ratio, rx2, so as to uniquely determine the other down- 
stream quantities u X 2, 0J pe 2, 0J C e2, and vt e 2 from the shock 
jump conditions for a magnetized two-fluid plasma consist- 
ing of electrons and ions with the equal bulk velocity and 
the equal number density. We adopt rx2 = 4.0 so that the 
thermal velocities of both downstream electrons and ions be- 
come much slower than the light speed. However, note that 
we can choose an arbitrary value for tt2- As can be seen 
later, the final cyclic reformation state does not depend on the 
choice of tt2- Since we performed the present simulation on 
a personal computer, we used a reduced ion-to-electron mass 
ratio r m = 100 for computational efficiency. With these 
parameters, we obtain the initial downstream quantities as 
^pe2 = 1-95, u> ce 2 = —0.19, u x 2 = 1.05, and vt e 2 = 7.55. 

We used 4096 cells for the upstream region and 8 192 cells 
for the downstream region. The grid spacing and time step 
of the present simulation are Ax/X e i — 1.0 and uj pe \At = 
0.01. We used 128 pairs of electrons and ions per cell in 
the upstream region and 512 pairs of electrons and ions per 
cell in the downstream region. It should be noted that such 
a few number of particles per cell is not enough to suppress 
the enhanced thermal fluctuations of particle-in-cell codes. 
In the present simulation, however, the numerical noises due 
to random motions of individual particles are substantially 
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Fig. 1. The transverse magnetic field B y as a function of position and 
time. The time and position are normalized by o) c i2 and pi = u x i/w c i2, 
respectively. The magnitude of magnetic field is normalized by Bo y i- 



3. Simulation Result 

Figure 1 shows the transverse magnetic field B y as a func- 
tion of position and time. The discontinuity exists at x = 
initially. However, the shock front shifts downstream due to 
the cyclotron motion of upstream ions that penetrate into the 
downstream region. A new shock front appears at x/pi ~ 1. 
The period from cj^i = to 7 is the transition epoch from 
the initial to the cyclic reformation state (see below). It 
should also be noted that the initial discontinuity causes the 
non-physical generation of electromagnetic waves at the on- 

OB 

set because Ampere's law = [Iq3 z is not satisfied at the 
initial state. However, the electromagnetic waves are per- 
fectly absorbed at both boundaries, and they never affect the 



T. UMEDA AND R. YAMAZAKI: A SHOCK-REST-FRAME MODEL 



3 




Fig. 2. Spatial profiles of the magnetic field energy density \fj,oB y \ 2 , ion-to-electron temperature ratio rx, number densities rt, thermal energy density 
components parallel and perpendicular to the ambient magnetic field TEu and TEj_, the drift energy density DE, and perpendicular kinetic energy 
versus position phase-space diagrams for electrons and ions at u) ci 2t = 38.1. The number densities are normalized by the initial upstream electron 
density n e i, and the energy densities are normalized by the initial upstream thermal energy density of electrons n e im e v^ el . The dashed lines show 
the downstream values obtained from the shock jump conditions for anisotropic plasmas with the downstream quantities v te \\ 2 = Vtel, r T\\2 = r Ti, 
rr±2 = 10, and the shock speed U /u x i = 0.0282. The perpendicular kinetic energies for electrons and ions are normalized by their upstream bulk 
energies, i.e., ^m e u^ 1 and ^"litt^ respectively. 



The typical processes for the self-reformation of perpen- 
dicular collisionless shocks are as follows (Quest, 1985; 
Lembege and Savoini, 1992; Schmitz et al., 2002a, b; 
Nishimura et al, 2003; Hada et al, 2003). A portion of 
ions are reflected back from the shock ramp. The shock 
foot region broadens toward the upstream region as the re- 
flected ions penetrate. However, the reflected ions return to 
the downstream region because of the cyclotron motion. The 
reflected ions can interact with upstream ions resulting in the 
self-formation of a new shock ramp in the old shock foot re- 
gion. Finally, a new shock front appears at x/pi ~ 1, and 
ions are again reflected toward the upstream region from the 
new shock ramp. 

The present simulation result shows that the timescale of 
the self reformation is almost equal to the cyclotron period 
of the downstream ions {u) C i2~t — 27r), which is in agree- 
ment with the previous studies (Quest, 1985, Lembege and 
Savoini, 1992; Schmitz et al, 2002a, b; Nishimura et al, 
2003; Scholer et al, 2003; Lee et al, 2004). Scholer et al 

(2003) showed that the self reformation process is not a com- 
putational artifact with the real ion-to-electron mass ratio, 
while the ratio of electron plasma-to-cyclotron frequency is 
smaller (oj pe i/ui ce i < 10). On the other hand, Lee et al. 

(2004) generated the self reformation of more than six cy- 
cles, while their mass ratio is much smaller (r m = 20). In the 
present shock rest frame, we have also confirmed the shock 
reformation process up to six cycles with r m = 100, but with 
less grid cells. 

In order to analyze the downstream condition in the cyclic 
reformation state, we plot in Figure 2 spatial profiles of the 



perature ratio perpendicular to the magnetic field rr±, the 
number density n, parallel and perpendicular thermal energy 
density components TE\\ and TE±, and the drift energy 
density DE for electron and ions at uj^t = 38.1. 

In the present shock-rest-frame model, a shock transition 
layer is self-consistently formed due to relaxation of the two 
plasmas with different quantities. The downstream region 
at the resulting cyclic reformation state is quite different 
from that of the initial state. For both electrons and ions, 
the spatial profiles of the thermal energy density component 
parallel to the ambient magnetic field are similar to those of 
number densities. This means that the downstream parallel 
temperatures for both electrons and ions become almost the 
same as those in the upstream region, i.e., T e ||i ~ T e || 2 — 
~ Ti\\i- On the other hand, the ratio of the ion-to- 
electron perpendicular temperature ratio is very large (~ 
10 2 ) in the transition region, and is typically tti_2 ~ 10 
in the downstream region (rx±2 = 5.0 ~ 20.0). 

As seen in Figure 1, the excited shock wave propagates 
slowly downstream with the roughly estimated shock speed 
U ju x \ = 0.0282. The downstream electron and ion number 
densities in the cyclic reformation state become smaller than 
those at the initial state, whereas the downstream electron 
and ion bulk velocities become faster than those at the initial 
state. We performed several additional runs with different 
system sizes and found that the simulation results — i.e., the 
shock speeds and the spatial profile of all physical quantities 
at an arbitrary time — are almost unchanged. In other words, 
the present system size is long enough to discuss the kinetic 
processes in the shock transition region. 
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the physical quantities fluctuate and are not spatially uni- 
form. We solved the shock jump conditions for anisotropic 
plasmas (Hudson, 1970) as a reference. Taking into ac- 
count the shock speed and the typical downstream quan- 
tities v te \\ 2 = v te i, r T \\ 2 = r T1 , and r T±2 = 10, the 
other downstream quantities are obtained as Lo pe2 = 1.71, 
w ce 2 = —0.146, u x2 = 1.33, and v te ±2 = 5.42. These 
quantities are plotted in Figure 2 with dashed lines. A differ- 
ence between the quasi-steady state by the shock jump con- 
ditions and the simulated downstream state is because of the 
dynamical shock reformation process. This might be another 
reason why we did not obtain the rigorous shock rest frame. 

For a given set of upstream parameters, the fluid shock 
jump conditions cannot give the downstream state uniquely. 
In the electron-ion fluid, the total plasma pressure is defined 
as the sum of the electron and ion pressures, P = n(T e +Ti), 
and the shock jump conditions allow us to take an arbitrary 
downstream temperature ratio, tt2- However, the value of 
rT2 is determined by the kinetic dynamics as seen in the 
simulation result. We performed several additional runs with 
different initial downstream values and confirmed that the 
shock speed in our reference frame and downstream thermal 
properties at the cyclic reformation state do not depend on 
the initial value of tt2 but on the other initial downstream 
parameters, such as the magnetic field, the number density, 
the bulk velocity, and the total plasma pressure. 

In the bottom panels of Figure 2 we show perpendicu- 
lar kinetic energy versus position diagrams for electrons and 
ions. We found that there exists a supra-thermal component 
of electrons at two local areas. The previous works reported 
the electron surfing acceleration with electrostatic solitary 
waves in the transition region (Hoshino and Shimada, 2002; 
Schmitz et al, 2002a, b). In the present simulation, we con- 
firmed the existence of solitary waves in the shock foot re- 
gion, where the maximum kinetic energy is about 30-fold 
more than the initial kinetic energy. On the other hand, an- 
other supra-thermal component of electrons due to the ion 
deceleration is seen at the overshoot, where the maximum ki- 
netic energy is much more (> 60KE e i) than that of the non- 
thermal electrons via the surfing mechanism. Although this 
process was also found in the previous simulations (Hoshino 
and Shimada, 2002; Schmitz et al, 2002a, b), the maximum 
kinetic energy in the present simulation is more than that in 
their simulations. We expect that this process becomes more 
significant with a larger mass ratio. 

4. Conclusion 

We have developed a shock-rest-frame model for full par- 
ticle simulations of perpendicular collisionless shocks based 
on the relaxation method used in the previous hybrid sim- 
ulations (e.g., Leroy et al, 1981, 1982). We reconfirmed 
both the formation of microscopic solitary structures due 
to the current-driven instability and cyclic reformation for 
a long time with a much smaller simulation domain. The 
shock-rest-frame model allows us to perform multidimen- 
sional full particle simulations of planar shocks more readily 
with current supercomputers. We are extending the present 
shock-rest-frame model to oblique shocks. Preliminary re- 
sults show that the present rest-frame model is very useful 
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